Shape  functions  for  three-dimensional  control- volume  mixed 
finite-element  methods  on  irregular  grids 

Richard  L.  Naff  a,  Thomas  F.  Russell  b  *  and  John  D.  Wilson  b  t 

aU.S.  Geological  Survey,  Denver,  CO  USA 

bUniversity  of  Colorado  at  Denver,  Denver,  CO  USA 

Numerical  methods  based  on  unstructured  grids,  with  irregular  cells,  frequently  require 
discrete  shape  functions  to  approximate  the  distribution  of  quantities  across  cells.  For 
control-volume  mixed  finite-element  methods,  vector  shape  functions  are  used  to  approx¬ 
imate  the  distribution  of  velocities  across  cells.  Previous,  two-dimensional  developments 
used  linear  shape  functions  to  interpolate  velocities  within  a  quadrilateral  cell.  For  ir¬ 
regular  hexahedral  cells  in  three  dimensions,  it  can  be  shown  that  linear  shape  functions 
cannot  exactly  represent  the  flux  distribution  across  a  cell  under  uniform  flow  conditions. 
As  a  result,  uniform  flow  cannot  be  exactly  simulated.  A  new  vector  shape  function  is 
proposed  for  use  with  irregular  hexahedral  cells  that  should  provide  for  a  more  accu¬ 
rate  velocity  approximation  within  a  cell.  This  velocity  shape  function  is  a  non-linear 
interpolator,  containing  quadratic  terms. 

1.  INTRODUCTION 

For  simulation  of  two-dimensional  flow  in  heterogeneous  porous  media,  it  has  been 
shown  that  mixed  methods,  and  in  particular  the  control-volume  mixed  finite-element 
(CVMFE)  methods,  are  often  the  most  accurate  methods  for  solving  for  the  velocity  field 
[1,2].  In  this  paper,  we  report  on  a  three-dimensional  velocity  shape  functions,  based 
on  covariant  vectors  for  a  mapping  to  a  unit  cube  and  for  use  with  irregular  hexahedral 
cells  [3].  The  three-dimensional  algorithm  described  herein  is  based  on  the  CVMFE 
methodology  as  developed  by  Cai  et  al.  [2]  for  the  simulation  of  Darcian  flow  in  two 
dimensions.  In  the  CVMFE  method,  the  domain  is  discretized  into  hexahedral  cells  that 
can  have  irregular  shapes,  allowing  for  the  modeling  of  complex  hydrogeological  systems. 
Shape  functions  serve  as  vector  basis  functions  to  interpolate  the  velocity  over  the  cell 
interiors.  Vector  test  functions  are  used  as  weighting  factors  for  integrating  the  Darcy 
relation  over  control  volumes  associated  with  cell  faces;  this  usage  can  be  viewed  as  an 
error  minimization  step  in  the  control  volume  technique  [4],  When  used  to  approximate 
the  Darcy  relation,  the  CVMFE  method  results  in  sets  of  discrete  equations  from  which 
bulk  fluxes  at  cell  faces  and  pressures  at  cell  centers  are  solved  for.  Shape  functions 
have  the  role  of  using  the  estimated  bulk  fluxes  at  the  cell  surface  to  approximate  the 
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Figure  1.  Left:  reference  cube  Q  with  edges  of  unit  length.  Right:  arbitrary  cell  Q 
from  discretization  with  vertex  locations  v0oo,  v10o,  v0i0,  v0oi,  v110,  v10i,  v011  and  vm 
indicated. 


velocity  in  the  cell  interior.  If  this  approximation  is  poor,  then  the  solution  obtained 
by  the  CVMFE  method  degrades.  The  velocity  shape  functions  for  three-dimensional 
logically  regular  meshes  proposed  in  [3]  should,  in  most  cases,  provide  a  reasonable  cell- 
velocity  estimate.  These  shape  functions  are  based  on  a  second-order  approximation  of 
flux  conservation  across  the  cell,  which  in  general  provides  for  a  non-linear  interpolation 
of  fluxes  and  the  velocities.  For  three-dimensional  flow  simulation  on  irregular  meshes, 
we  believe  that  this  function  offers  advantages  over  the  linear  shape  function  currently 
used  in  two-dimensional  simulations  [2,5]. 

2.  BASIC  EQUATIONS 

The  numerical  method  outlined  in  this  paper  is  based  on  the  following  steady  flow 
equations  applied  within  a  three-dimensional  domain  fh 

V-q  =  W(x,y,z)  and  q  =  -K  (x,y,z)Vp/p,  (x,y,z)<EVL.  (1) 

Here  q  is  the  specific  discharge  vector,  W ( x ,  y.  z )  is  a  source  term,  p  is  the  pressure, 
K (x,y,z)  is  the  hydraulic  conductivity  tensor  and  p  is  a  general  viscosity  term.  On  the 
surface  of  the  domain,  boundary  conditions  can  consist  of  specified  fluxes  over  dQ j  or 
specified  pressures  over  dQp.  For  the  CVMFE  method,  the  hydraulic  conductivity  tensor 
is  inverted,  causing  the  Darcy  relation  in  (1)  to  be  written  as 

Vp  =  — pK^q.  (2) 

The  mixed-method  development  herein  uses  the  continuity  relation  from  (1)  and  the 
inverted  Darcy  relation  from  (2)  as  a  basis  of  the  numerical  approximation. 

The  domain  is  discretized  using  a  logically  regular  mesh  of  hexahedral  cells,  as  defined 
subsequently.  This  mesh  may  be  irregular  in  construction  to  the  point  where  the  primary, 
bounding  faces  of  cells  are  not  planar.  Any  given  hexahedral  cell  Q  can  be  described 
from  the  location  of  its  vertices  at  v0oo,  v10o,  v0io,  v0oi,  v110,  v10i,  v0n  and  vm,  where 
va/37  =  (xa$fe  yafir>  za/3-/)  (Figure  1).  A  trilinear  mapping  exists  such  that  the  hexahedral 
cell  Q  is  the  image  of  a  regular  reference  hexahedron,  Q,  consisting  of  a  unit  cube  with 
fixed  vertices  at  (0,  0,  0),  (1,  0,  0),  (0, 1,  0),  (0,  0, 1),  (1, 1,  0),  (1,  0, 1),  (0, 1, 1)  and  (1, 1, 1), 


on  a  point  by  point  basis  (see  Figure  1).  This  mapping  allows  that,  for  any  reference 
location  r  =  ( x ,  y.  5)  within  Q,  the  equivalent  location  r  =  ( x ,  y,  z )  within  Q  is  obtained 
from  the  following  expression: 

r  =  v0  +  vaf  +  vby  +  vcz  +  vdxy  +  vexz  +  vfyz  +  vgxyz  (3) 

where  v0  =  v000,  va  =  v100  -  v0,  =  v0i0  -  v0,  vc  =  v00i  -  v0,  vd  =  vno  -  v0  -  va  -  v6, 

Ve  =  Vioi-V0-Va-Vc,  V/  =  V0n-Vo-Vft-Vc,  Vfl  =  Vm  V0  Va  ~Vf)  vc  Ve  Vj. 

(Also  see  [2,5].)  An  irregular  mesh  of  hexahedra  Q  has  a  reference  mesh  of  regular  cubes; 
such  meshes  are  referred  to  as  logically  regular.  Directions  associated  with  the  reference 
mesh  are  referred  to  as  the  logical  x.  y  and  z  directions.  Note  that,  should  x  be  fixed  in 
f,  then  a  face  normal  to  the  x  direction  within  Q  is  determined.  This  face  is  mapped, 
via  (3),  into  an  equivalent  face  within  Q.  Because  0  <  x  <  1,  such  faces  are  referred 
to  as  intermediate  interior  cell  faces.  Covariant  vectors,  defined  as  X($,  5)  =  dr/dx, 
Y (x,z)  =  dr/dy  and  Z (x,y)  =  dr/dz,  allow  for  definition  of  the  geometry  of  Q.  The 
volumetric  Jacobian  J  for  passing  from  Q  to  Q  is  simply  [6] 

J(x,y,z)  =  X(y,z)  ■  (Y(x,z)  x  Z(x,y))  (4) 

Surface  Jacobians  are  similarly  defined  for  faces  in  the  logical  x ,  y  and  z  directions  [3]. 

3.  VELOCITY  SHAPE  FUNCTIONS 

The  shape  functions  proposed  in  [3]  offer  advantages  for  general  hexahedral  cells;  we 
present  them  here.  The  CVMFE  method  uses  velocity  basis  vectors,  or  shape  functions, 
to  approximate  the  Darcy  velocity  q  in  an  arbitrary  cell  Q  from  bulk  fluxes  at  cell  faces. 
Note  that  “flux”  is  defined  herein  to  be  the  total  volumetric  discharge  through  a  cell  face, 
as  opposed  to  the  volumetric  discharge  per  unit  area.  The  fluxes  at  every  cell  face  are 
denoted  as  fx 0,  fx i,  fy o,  fy i,  fz o  and  fzi  where,  for  instance,  fx 0  is  the  flux  through  the 
face  at  x  =  0.  The  shape  functions  for  the  cell  faces  are  denoted  as  vx 0,  i/xi,  vy0,  isyl,  vz0 
and  vzX.  Allowing  Vc  to  be  this  approximation  of  q,  then  its  cell- wise  representation  is 

Vc  =  fxO^xO  +  fx \vx\  +  fyOVyO  +  fylVyi  +  fz  0^20  +  f zlv  zl-  (5) 

The  Darcy  velocity  q  is  represented  by  such  a  relation  for  every  cell  in  the  mesh.  For 
general  hexahedral  cells,  the  shape  functions  are  assumed  to  have  the  following  form 
[3]: 

*Vo  =  j[(l  -  x)X/3x0  -  x{\  -  x)X/3x2rxx „  -  y(  1  -  y)Y(3y2ryx0  - 1(  1  -  z)Z/3z2rzx0\,  (6a) 
vx i  =  -j\xX/3xl  -  x{\  -  x)Xfix2rxx i  -  y(  1  -  y)Y (3y2ryxl  -  z(  1  -  z)Z(3z2rzxl\,  (6b) 

vyo  =  j[(l  -  y)Y(3y0  -  x(l  -  x)X(3x2rxy 0  -  y{  1  -  y)Y /3y2ryy0  -  z(l  -  z)Zf3z2rzy0\,  (6c) 
vyi  =  j[y Y0yl  -  x(l  -  x)X/3x2rxyl  -  y(  1  -  y)Y f3y2ryyl  -  5(1  -  z)Z(3z2rzyl\,  (6d) 

vz0  =  j[(l  -  z)Z/3z0  -  x(l  -  x)X/3x2rxz 0  -  y(  1  -  y)Y /3y2ryz0  -  5(1  -  z)Z/3z2rzz0\,  (6e) 

vz i  =  j[zZ/3zl  -  x(l  -  x)Xj3x2rxzi  -  y(  1  -  y)Y fiy2ryzl  -  5(1  -  z)Z/3z2rzzl],  (6f) 


where,  for  k  =  x,  y,  or  z,  and  ^  =  0,1,  the  r  factors  in  (6a)  -  (6/)  become 
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Here,  J(x,y,z)  is  the  volumetric  Jacobian  (4).  The  factors 
are  given  by 
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where  Ax(x ),  Ay(y)  and  Az(z)  are  surface  areas  of  intermediate  interior  faces: 
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Similarly,  the  expressions  for  /3x2,  /3J/2  and  /Jz2  are 
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Just  as  the  cross  products  in  (8)  define  cell  surfaces,  the  cross  products  contained  in 
define  a  set  of  secondary  surfaces.  The  vectors  W2x,  W 2y  and  W2z  take  the  form 

W2a,  =  Ya*  (1/2)  x  Za*  (1/2) ,  W 2y  =  Z2j  (1/2)  x  X2j  (1/2) , 
w2z  =  X2z  (1/2)  x  Y2z  (1/2) . 


(7) 

0,1 
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The  vectors  associated  with  secondary  surfaces  are:  X2z(y)  =  X|^=1  —  X|^=0,  X2j/(i)  = 
X|  jj=i  X|  jj=o,  Y2x(z)  Y  |^=i  Y  | £=0 5  Yfiz (d )  Y 1 2=  i  Y  \z=o ,  Z2a,(y)  Z  x=  i  Z  \x=q 
and  Z 2y(x)  =  Z|^=1  —  Z|^=0.  It  is  shown  in  [3]  that  (6a)-(6/)  are  second-order  correct  in 
approximation.  Note  that  if  shape  functions  were  approximated  with  linear  forms  then 
only  the  first  term  in  each  of  (6a)-(6/)  would  be  present  [4],  If  the  mesh  elements  are 
parallelepipeds,  then  the  quadratic  terms  in  (6a)-(6/)  will  all  be  null;  for  such  regular 
meshes,  there  is  no  advantage  to  these  shape  functions. 

Using  (5)  to  represent  q  in  (1),  the  following  discrete  continuity  equation  results  when 
both  sides  of  the  continuity  equation  in  (1)  are  integrated  over  Q: 

fx  1  -  fxo  +  fyl  -  fyo  +  fzi  ~  fzo  =  W (x,  y,  z)  dxdydz.  (12) 

J  Q 

This  result  follows  from  Gauss’  divergence  theorem  and  properties  of  shape  functions  [3]. 
A  set  of  discrete  Darcy  relations  are  obtained  by  integrating  (2)  against  a  test  function 
over  subdomains  spanning  all  cell  faces  of  the  discretized  domain  [2-4];  an  equation  similar 
to  (5)  is  used  to  approximate  q  in  these  subdomains.  The  discrete  Darcy  relation  and  the 
discrete  continuity  equation  (12),  applied  to  every  cell  within  the  domain,  form  the  sets 
of  equations  which  are  the  basis  of  the  CVMFE  method. 
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Figure  2.  a.  Hexahedral  cell  in  form  of  a  truncated  pyramid.  Uniform  flow  parallel  to 
vertical  z  axis,  as  indicated  by  arrow,  is  assumed.  s0  and  si  indicate  the  length  of  the 
lower  and  upper  edges,  respectively,  b.  Hexahedral  cell  in  form  of  a  tent  with  curved 
roof.  Uniform  flow  parallel  to  horizontal  x  axis,  as  indicated  by  arrow,  is  assumed.  Dash- 
dot  lines  show  maxima  and  minima  of  saddle  surface  forming  cell  face  2  =  1.  For  both 
figures,  vertices  v0oo  and  vm  are  indicated;  see  Figure  1  for  orientation. 

4.  UNIFORM  FLOW  TESTS 

For  general  hexahedral  cells  and  under  uniform  flow,  linear  shape  functions  of  the 
Piola  type  will  not  give  the  proper  flux  through  an  arbitrary  cell  face.  As  a  simple 
example,  consider  the  “truncated  pyramid”  hexahedron  Q  (Figure  2a),  whose  bottom 
face  at  z  =  (s0  —  S\)/2  ( 2  =  0)  is  a  square  of  side  s0  given  by  — s0/2  <  x  <  s0/ 2, 
—s0/ 2  <  y  <  s0/2  (with  area  Sq),  and  whose  top  face  at  z  =  (si  —  s0)/2  ( 2  =  1)  is  the 
square  of  side  si  given  by  —  S\/2  <  x  <  Si/2,  —  S\/2  <  y  <  S\/2  (with  area  sf): 

x  =  a(2x  —  l)/2,  y  =  a(2y  —  l)/2,  z  =  (2  —  l/2)(si  —  s0)  (13) 

where  a  =  s0  +  2(si  —  s0).  The  four  lateral  faces  are  planar  trapezoids  (see  Figure  2a). 
For  a  unit  vertical  uniform  flow  q  =  (0,  0, 1),  the  exact  flux  across  horizontal  interior  faces 
for  fixed  2  is  the  cross-sectional  area, 

/*(2)  =  [(1  -  2)s0 +  2si]2.  (14) 

A  linear  shape  function,  while  producing  the  correct  fluxes  at  the  top  and  bottom,  will 
approximate  the  flux  at  a  face  2  as 

W)  ~  (1  -  £)/2(0) +  2/2(l)  =  (1-2)sq  +  2s2.  (15) 

A  quadratic  flux  interpolation,  of  the  type  contained  in  (6a)-(6/),  can  deal  exactly  with 
this  example.  For  uniform  flow  in  general,  (5)  with  these  shape  functions  gives  an  exact 
cell  velocity  interpolation  if  the  primary  cell  faces  and  secondary  surfaces  (defined  after 
(10))  are  all  planar  [3].  Mesh  construction  can  usually  produce  planar  primary  faces,  but 


Figure  3.  L2  norm  for  mesh  composed  of 
cells  with  quasi-random  planar  faces;  see 
text  for  h.  Solid  line:  non-linear  shape 
function;  Dashed  line:  linear  shape  func¬ 
tion. 


Figure  4.  L2  norm  for  mesh  composed 
of  cells  with  random  vertices;  see  text  for 
h.  Solid  line:  non-linear  shape  function; 
Dashed  line:  linear  shape  function. 


planar  secondary  surfaces  can  be  more  problematic;  in  particular,  secondary  surfaces  can 
be  non-planar  even  when  all  exterior  cell  faces  are  planar.  However,  our  experience  has 
been  that,  even  if  the  secondary  surfaces  are  non-planar,  quite  reasonable  results  can  be 
obtained  as  long  as  the  primary  faces  are  planar  [3].  Planar  secondary  surfaces  can  be 
forced  by  insisting  that  each  cell  have  at  least  one  pair  of  parallel  opposite  faces. 

Use  of  shape  functions  (6 a)-(Qf)  becomes  problematic  when  the  primary  faces  are  non- 
planar;  an  example  of  a  non-planar  primary  face  is  the  “tent  with  curved  roof”  (Figure 
2b)  with  unit  square  base  at  z  =  0  (z  =  0),  vertical  lateral  faces,  and  roof  height  h0  at 
two  opposite  vertices  (0,  0,  h0),  (1, 1,  h0 )  and  height  hi  at  the  other  two  vertices: 

x  =  x,  y  =  y.  z  =  (h0[(l  -  x)(l  -  y)  +  xy\  +  hi[(l  -  x)y  +  x (1  -  y)])z.  (16) 

In  (6o)-(6/)  the  covariant  vectors  X,  Y  and  Z  provide  the  directional  underpinnings  for 
the  shape  functions.  However,  at  face  z  =  1,  X,  Y,  and  (1  —  z) Z  all  have  null  normal 
component,  and  fz i  =  0  in  (5)  for  q  =  (1,  0,  0).  Thus  Vc  in  (5)  has  null  normal  component 
at  all  points,  while  q  does  not.  Hence  (6a)-(6/)  cannot  interpolate  uniform  q  exactly. 

Two  sets  of  simulations  were  performed  to  demonstrate  these  effects;  in  all  cases,  meshes 
with  8x8x8  cells  were  used.  In  the  first,  irregular  cells  with  planar,  quasi-random  faces 
were  generated,  starting  with  random  vertices  on  three  adjoining  exterior  faces  of  O. 
These  vertices  were  created  from  a  regular,  two-dimensional  mesh  on  each  face  by  adding 
to  each  regular  vertex  location  a  uniform  random  variable  from  the  interval  (—h/ 2,  h/2). 
Starting  at  the  corner  where  the  three  exterior  faces  meet,  cells  were  created  by  connecting 
vertices  of  three  existing  cell  faces  to  an  interior  vertex,  requiring  that  the  new  cell  faces 
forming  this  connection  be  planar.  This  results  in  interior  cells  with  planar  faces,  none 


Figure  5.  Example  of  three  exterior 
faces  of  a  4x4x4  domain  that  share  vertex 
(0,0,0).  Fold  up  the  left  ( yz )  face  along  the 
y-axis  and  the  bottom  (xz)  face  along  the 
x-axis.  Randomness  of  vertex  locations  is 
exaggerated. 


Figure  6.  L-2  norm  for  mesh  composed  of 
cells  with  quasi-random  planar  faces;  see 
text  for  (As).  Solid  line:  non-linear  shape 
function;  Dashed  line:  linear  shape  func¬ 
tion;  Circle:  d  is  one  domain  length;  Cross: 
d  is  two  domain  lengths. 


of  which  are  parallel.  The  secondary  surfaces  of  cells  in  this  grid  are  rarely,  if  ever, 
planar.  Error  in  simulating  uniform  flow  parallel  to  the  x  direction  with  meshes  of  this 
type  was  quantified  using  the  L-2  norm  [4],  Figure  3  is  a  depiction  of  the  L-2  norm  results 
for  this  case;  increasing  h  caused  successive  meshes  to  be  increasingly  irregular.  These 
results  demonstrate  the  superiority  of  the  non-linear  shape  function  over  the  linear  shape 
function  for  simulating  uniform  flow.  As  the  secondary  surfaces  associated  with  vectors 
X2j/,  X2„  Y2x  ,  Y2z,  Z2x  and  Z2y  are  non-planar,  the  non-linear  shape  function  provides 
an  improved  result,  but  not  an  exact  result.  In  the  second  set  of  simulations  (Figure 
4),  the  mesh  was  created  by  perturbing  all  the  vertex  locations  of  a  regular  mesh  with 
a  uniform  random  variable;  the  vertices  on  dQ  were  perturbed  only  within  the  plane 
forming  the  surface.  This  procedure  results  in  meshes  where  none  of  the  cell  faces  interior 
to  are  planar;  uniform  flow  in  the  x  direction  was  simulated.  The  results  shown  in 
Figure  4  indicate  no  particular  advantage  for  the  non-linear  shape  function.  As  this  set  of 
simulations  primarily  tests  the  non-planar  aspect  of  the  primary  cell  faces,  these  results 
confirm  the  observations  regarding  the  ‘tent  with  curved  roof”  cell  depicted  in  Figure  2b. 

5.  NON-UNIFORM  FLOW  TESTS 

For  non-uniform  flow,  a  dipole  test  was  devised  whereby  a  cubic  domain  (four  units  on 
a  side)  was  located  midway  between  a  source  and  a  sink  of  equal  strength;  the  domain 
was  oriented  such  that  the  edge  y  =  0,  z  =  0  coincided  with  a  straight  line  between 
the  poles.  The  distance  from  a  dipole  to  the  closest  domain  corner  d  was  allowed  to  be 
either  one  or  two  domain  lengths.  Flux  boundary  conditions  for  the  simulations  were 


obtained  by  integrating  the  normal  component  of  the  dipole  velocities  over  the  exterior  of 
the  domain;  these  boundary  flux  estimates  were  generally  very  accurate.  Quasi-random 
meshes  with  planar  faces  were  generated  by  the  same  method  as  before,  except  that 
now  all  cell  faces  perpendicular  to  the  x  direction  are  parallel  (Figure  5);  vertices  on  the 
exterior  faces  were  randomized  in  the  y  and  z  directions  with  a  uniform  random  variable 
from  the  interval  (—0.1,  0.1).  This  procedure  results  in  planar  secondary  surfaces  for  all 
cells  [3].  Meshes  with  4x4x4,  8x8x8,  16x16x16  and  32x32x32  cells  were  generated 
in  this  manner,  corresponding  to  mean  discretization  lengths  (As)  of  1,  1/2,  1/4  and  1/8 
respectively.  The  results  from  these  simulations,  characterized  by  their  L2  norms,  are 
depicted  in  Figure  6.  These  results  indicate  that,  at  the  coarsest  discretization,  the  error 
from  approximating  non-uniform  flow  overwhelms  any  corrective  effect  of  the  non-linear 
shape  function  on  the  simulation.  This  error  results  because  the  flux  estimates  on  the  cell 
surfaces  produced  by  CVMFE  are  average  or  bulk  estimates.  However,  with  refinement 
of  the  grid  (decreasing  (As)),  the  non-linear  shape  function  outperforms  the  linear  shape 
function;  with  (As)  =  1/8,  an  order  of  magnitude  difference  can  be  observed.  This 
results,  in  large  part,  because  the  linear  shape  functions  cease  to  provide  improvement  in 
performance  with  increased  grid  refinement  when  the  discretization  is  quasi-random  (see 
[3]  for  a  more  theoretical  discussion  of  the  linear /non-linear  shape  function  behavior). 

6.  CONCLUDING  REMARKS 

In  [3]  we  note  that  the  convergence  of  CVMFE,  with  shape  functions  (6a)  (6/),  is 
second  order.  Here  it  is  demonstrated  that  the  non-linear  shape  functions  can  provide  for 
better  resolution  of  uniform  flow  with  general  meshes,  relative  to  linear  shape  functions, 
provided  that  primary  cell  faces  are  planar.  For  non-uniform  flow,  the  level  of  mesh 
refinement  is  all  important;  only  when  a  general  grid  is  relatively  fine  will  the  non-linear 
aspects  of  the  proposed  shape  functions  contribute  significantly. 
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